# Who Do You Trust? - Daniel Goldstein and Johannes Wiedemann
# R-script for Analysis (Web Appendix)


#====================================================================================================


# load packages

library(tidyverse)
library(readxl)
library(lfe)
library(interactions)
library(estimatr)
library(openxlsx)
library(texreg)
library(dotwhisker)
library(stargazer)
library(cobalt)
library(stringi)



#=======================================================================
# read in data

# Change the working directory as appropriate
setwd("~/../GoldsteinWiedemann_supplementary")





load(file="Data and Data Preparation/GW_2020_MainDataSet_122220.RData") 
load(file="Data and Data Preparation/GW_2020_DataSetDaily2019_122220.RData") 
load(file="Data and Data Preparation/GW_2020_MRPDataSet_122220.RData")




#=====================================================================================================
#=====================================================================================================





#====================================================================================================
# Appendix-Analysis starts here




#====================================================================================================
# TABLES (Appendix)


#=====================================================================================================
# Create Summary Table

data_weekly_series_1 %>% ungroup()%>% mutate(PerCapitaInclog=log(PerCapitaInc),PopDensity2010log=log(PopDensity2010+1)) %>% select(cmi_weekly,cmi_weekly_2019,shelter_2020,sheltered_in_place_2019,gop_per_2016,sk2017_alt,Turnout12and16,Republican_gov_2020,analysis_value,PerCapitaInclog,UnempRate2018,Ed5CollegePlusPct,WhiteNonHispanicPct2010,BlackNonHispanicPct2010,AsianNonHispanicPct2010,Age65AndOlderPct2010,PopDensity2010log,RuralUrbanContinuumCode2013,PctEmpAgriculture,PctEmpManufacturing,PctEmpServices) %>%
  as.data.frame()%>%
  stargazer(summary = TRUE,table.placement = "h",float = FALSE,header=FALSE,median = TRUE,style = "aer",font.size="small",covariate.labels= c("CMI 2020","CMI 2019","Pct Shelter 2020","Pct Shelter 2019","Republican Per.","Social Capital","Turnout","Repub. Governor","Vaccination Per.","Income p.c.","Unemployment","Education","Pct White","Pct Black","Pct Asian","Pct over 65","Pop. Dens.","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv."))




#============================================================================


# Create Regression Tables for Stay-at-Home Effects by Turnout


panel_9 <- felm(cmi_weekly~post_treatment3*Turnout12and16|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_b <- felm(cmi_weekly~post_treatment3*Turnout12and16+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_c <- felm(cmi_weekly~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_d <- felm(cmi_weekly~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_e <- felm(cmi_weekly~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_f <- felm(cmi_weekly~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_9,panel_9_b,panel_9_c,panel_9_d,panel_9_e,panel_9_f),stars = c(0.01,0.05,0.1),table = FALSE,custom.coef.names = c("Order","Turnout per.","Order x Turnout per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "Mobility after Shelter Order by Turnout (Continuous Measure)",caption.above = TRUE,lyx = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,use.packages = FALSE)


#====================================================================================================




# Create Regression Tables for Baseline Treatment Effects, with Pct Shelter instead of CMI as Dependent Variable


panel_1 <- felm(shelter_2020~post_treatment3|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_b <- felm(shelter_2020~post_treatment3+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_c <- felm(shelter_2020~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_d <- felm(shelter_2020~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_e <- felm(shelter_2020~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_f <- felm(shelter_2020~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_1,panel_1_b,panel_1_c,panel_1_d,panel_1_e,panel_1_f),stars = c(0.01,0.05,0.1),table = FALSE,custom.coef.names = c("Order","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Pct Sheltering after Shelter Orders",caption.above = TRUE,lyx = TRUE,scalebox = 0.7,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,use.packages = FALSE)



#====================================================================================================


# Create Regression Tables for Partisan Effects, with Pct Shelter instead of CMI as Dependent Variable

panel_3 <- felm(shelter_2020~post_treatment3*gop_per_2016|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_b <- felm(shelter_2020~post_treatment3*gop_per_2016+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_c <- felm(shelter_2020~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_d <- felm(shelter_2020~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_e <- felm(shelter_2020~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_f <- felm(shelter_2020~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_3,panel_3_b,panel_3_c,panel_3_d,panel_3_e,panel_3_f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Republican Per.","Order x Rep. Per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Sheltering after Shelter Order by Party (Continuous Measure)",caption.above = TRUE,lyx = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE)


#====================================================================================================

# Create Regression Tables for Governor Partisanship and County Partisanship Effects, with Pct Shelter instead of CMI as Dependent Variable

panel_14 <- felm(shelter_2020~post_treatment3*Republican_gov_2020*gop_per_2016|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_b <- felm(shelter_2020~post_treatment3*Republican_gov_2020*gop_per_2016+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_c <- felm(shelter_2020~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_d <- felm(shelter_2020~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_e <- felm(shelter_2020~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_f <- felm(shelter_2020~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_14,panel_14_b,panel_14_c,panel_14_d,panel_14_e,panel_14_f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Republican Governor","Republican per.","Order x Republican Gov.","Order x Rep. per.","Rep. Gov x Rep. per.","Order x Rep. Gov. x Rep. per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Sheltering after Shelter Order by Partisanship of Governor and County (Continuous Measure)",caption.above = TRUE,lyx = TRUE,scalebox = 0.4,include.rsquared=FALSE,custom.note = "",include.adjrs=FALSE,include.groups=FALSE)




#====================================================================================================

# Create Regression Tables for Social Capital Effects, with Pct Shelter instead of CMI as Dependent Variable

panel_5 <- felm(shelter_2020~post_treatment3*sk2017_alt|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_b <- felm(shelter_2020~post_treatment3*sk2017_alt+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_c <- felm(shelter_2020~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_d <- felm(shelter_2020~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_e <- felm(shelter_2020~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_f <- felm(shelter_2020~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_5,panel_5_b,panel_5_c,panel_5_d,panel_5_e,panel_5_f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Social Capital","Order x Soc. Cap.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Sheltering after Shelter Order by Social Capital (Continuous Measure)",caption.above = TRUE,lyx = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE)



#====================================================================================================

# Create Regression Tables for Social Capital and Partisanship Effects, with Pct Shelter instead of CMI as Dependent Variable

panel_15 <- felm(shelter_2020~post_treatment3*sk2017_alt*gop_per_2016|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_b <- felm(shelter_2020~post_treatment3*sk2017_alt*gop_per_2016+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_c <- felm(shelter_2020~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_d <- felm(shelter_2020~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_e <- felm(shelter_2020~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_f <- felm(shelter_2020~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_15,panel_15_b,panel_15_c,panel_15_d,panel_15_e,panel_15_f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Social Capital","Republican Per","Order x Soc. Cap.","Order x Rep.","Soc. Cap. x Rep.","Order x Soc. Cap. x Rep.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Sheltering after Shelter Order by Social Capital and County Partisanship",caption.above = TRUE,lyx = TRUE,scalebox = 0.4,include.rsquared=FALSE,custom.note = "",include.adjrs=FALSE,include.groups=FALSE)


#====================================================================================================

# Create Regression Tables for Vaccination Rate Effects, with Pct Shelter instead of CMI as Dependent Variable

panel_20 <- felm(shelter_2020~post_treatment3*analysis_value|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_b <- felm(shelter_2020~post_treatment3*analysis_value+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_c <- felm(shelter_2020~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_d <- felm(shelter_2020~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_e <- felm(shelter_2020~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_f <- felm(shelter_2020~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])



texreg(list(panel_20,panel_20_b,panel_20_c,panel_20_d,panel_20_e,panel_20_f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Vaccination Per.","Order x Vacc. per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Sheltering after Shelter Order by Vaccination Rate",caption.above = TRUE,lyx = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE)





#====================================================================================================

# Create Regression Tables for Turnout Effects, with Pct Shelter instead of CMI as Dependent Variable

panel_9 <- felm(shelter_2020~post_treatment3*Turnout12and16|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_b <- felm(shelter_2020~post_treatment3*Turnout12and16+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_c <- felm(shelter_2020~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_d <- felm(shelter_2020~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_e <- felm(shelter_2020~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+sheltered_in_place_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_9_f <- felm(shelter_2020~post_treatment3*Turnout12and16+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


texreg(list(panel_9,panel_9_b,panel_9_c,panel_9_d,panel_9_e,panel_9_f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Turnout per.","Order x Turnout per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","Shelter 2019"),caption = "Sheltering after Shelter Order by Turnout (Continuous Measure)",caption.above = TRUE,lyx = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE)






#================================================================


# Create Regression Table with Placebo Treatment Effects

data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac1=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-11),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac2=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-18),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac3=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(3),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac4=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(10),1,0,missing = 0))

panel_plac1 <- felm(cmi_weekly~post_treatment3|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac2 <- felm(cmi_weekly~post_treatment_plac1|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac3 <- felm(cmi_weekly~post_treatment_plac2|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac4 <- felm(cmi_weekly~post_treatment_plac3|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac5 <- felm(cmi_weekly~post_treatment_plac4|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)


texreg(list(panel_plac1,panel_plac2,panel_plac3,panel_plac4,panel_plac5),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Order Lag","Order Lag 2","Order Lead","Order Lead 2"),caption = "Mobility after Shelter Orders - Placebo Checks",caption.above = TRUE,lyx = TRUE,scalebox = 0.75,include.rsquared=FALSE,custom.note = "",include.groups=FALSE)




#================================================================



# Create Regression Table with Placebo Treatment Effects (with control variables)


data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac1=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-11),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac2=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-18),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac3=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(3),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac4=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(10),1,0,missing = 0))

panel_plac1 <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac2 <- felm(cmi_weekly~post_treatment_plac1+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac3 <- felm(cmi_weekly~post_treatment_plac2+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac4 <- felm(cmi_weekly~post_treatment_plac3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_plac5 <- felm(cmi_weekly~post_treatment_plac4+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)


texreg(list(panel_plac1,panel_plac2,panel_plac3,panel_plac4,panel_plac5),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,custom.coef.names = c("Order","Pct over 65","Pop. Dens.","Income p.c.","Unemployment","Education","Pct White","Pct Black","Pct Asian","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Order Lag","Order Lag 2","Order Lead","Order Lead 2"),caption = "Mobility after Shelter Orders - Placebo Checks",caption.above = TRUE,lyx = TRUE,scalebox = 0.75,include.rsquared=FALSE,custom.note = "",include.groups=FALSE)








#================================================================




# Sensitivity Analysis for Urban-Rural Differences



sens3a <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1[data_weekly_series_1$RuralUrbanContinuumCode2013>=7,]) # 1077 these are rural or small urban populations (<2,500); away from metro areas
sens3b <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1[data_weekly_series_1$RuralUrbanContinuumCode2013<7,]) # these are urban populations
sens3c <- felm(cmi_weekly~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1[data_weekly_series_1$RuralUrbanContinuumCode2013>=7,]) # 1077 these are rural or small urban populations (<2,500); away from metro areas
sens3d <- felm(cmi_weekly~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1[data_weekly_series_1$RuralUrbanContinuumCode2013<7,]) # these are urban populations
sens3e <- felm(cmi_weekly~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1[data_weekly_series_1$RuralUrbanContinuumCode2013>=7,]) # 1077 these are rural or small urban populations (<2,500); away from metro areas
sens3f <- felm(cmi_weekly~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1[data_weekly_series_1$RuralUrbanContinuumCode2013<7,]) # these are urban populations

texreg(custom.model.names = c("Rural","Urban","Rural","Urban","Rural","Urban"),list(sens3a,sens3b,sens3c,sens3d,sens3e,sens3f),stars = c(0.01,0.05,0.1),table = FALSE,use.packages = FALSE,caption = "Mobility after Shelter Orders - Placebo Checks",caption.above = TRUE,lyx = TRUE,scalebox = 0.75,include.rsquared=FALSE,custom.note = "",custom.coef.names = c("Order","Pct over 65","Pop. Dens.","Income p.c.","Unemployment","Education","Pct White","Pct Black","Pct Asian","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Republican Per.","Order x Rep. Per.","Social Capital","Order x Soc. Cap."),include.groups=FALSE)




#===========================================================================================================
#===========================================================================================================
# FIGURES (Appendix)


#===========================================================================================================


# Create sanity check on mobility data: can we detect changes in mobility around Thanksgiving 2019?

cuebiq_2019_daily %>% group_by(as.Date(ref_dt)) %>%mutate(mean_cmi=mean(cmi,na.rm=TRUE))%>% filter(as.Date(ref_dt)>"2019-11-22")%>%filter(as.Date(ref_dt)<"2019-12-01")%>%mutate(Date=as.factor(as.Date(ref_dt)))%>%
  mutate(Thanksgiving=ifelse(as.Date(ref_dt)=="2019-11-28",1,NA))%>%
  ggplot(aes(x=Date,y=cmi))+
  theme_bw()+
  labs(
    y = "CMI", x = "Date") + 
  geom_label(aes(label="Thanksgiving",x="2019-11-28",y=4.7),color="red")+
  geom_boxplot()

#ggsave("Figures PoP Appendix/CMIValidity_final.png")


#================================================================


# Create Parallel Trends plots for Mobility Trends by Treatment Week and States 



plot1 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_12_pre=ifelse(TreatmentWeek_num2==12&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_12_post=ifelse(TreatmentWeek_num2==12&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(early_late_never==2&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  mutate(State=as.factor(State_abb)) %>%
  ggplot(aes(x=week_num,color=State)) +
  stat_summary(aes(y=group_12_pre,linetype="solid",color=State),geom = 'line') +
  stat_summary(aes(y=group_12_post,linetype="dashed",color=State),geom='line') +
  stat_summary(aes(y=pre_3,linetype="longdash"),geom='line',color="black",lwd=1.4) +
  stat_summary(aes(y=pre_2,linetype="dotted"),geom='line',color="black",lwd=1.4) +
  geom_vline(aes(xintercept = 11),color="grey",linetype="dashed",lwd=1.3)+
  theme_minimal()+
  theme(legend.position = "bottom",legend.title = element_text(color = "black", size = 9),legend.text = element_text(color = "black", size = 7))+
  guides(size = FALSE)+
  xlim(8,16)+
  scale_color_discrete(name="")+
  scale_linetype_discrete(name="",labels=c("Never","Late"),breaks=c("longdash","dotted"))+
  xlab("Week Number")+
  ylab("Mobility")+
  ggtitle("Week 12")


plot2 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_13_pre=ifelse(TreatmentWeek_num2==13&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_13_post=ifelse(TreatmentWeek_num2==13&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(early_late_never==2&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  mutate(State=as.factor(State_abb)) %>%
  ggplot(aes(x=week_num,color=State)) +
  stat_summary(aes(y=group_13_pre,linetype="solid",color=State),geom = 'line') +
  stat_summary(aes(y=group_13_post,linetype="dashed",color=State),geom='line') +
  stat_summary(aes(y=pre_3,linetype="longdash"),geom='line',color="black",lwd=1.4) +
  stat_summary(aes(y=pre_2,linetype="dotted"),geom='line',color="black",lwd=1.4) +
  geom_vline(aes(xintercept = 12),color="grey",linetype="dashed",lwd=1.3)+
  theme_minimal()+
  theme(legend.position = "bottom",legend.title = element_text(color = "black", size = 9),legend.text = element_text(color = "black", size = 7))+
  guides(size = FALSE)+
  xlim(8,16)+
  scale_linetype_discrete(name="",labels=c("Never","Late"),breaks=c("longdash","dotted"))+
  scale_color_discrete(name="")+
  xlab("Week Number")+
  ylab("Mobility")+
  ggtitle("Week 13")

plot3 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_14_pre=ifelse(TreatmentWeek_num2==14&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_14_post=ifelse(TreatmentWeek_num2==14&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(TreatmentWeek_num2==15&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  mutate(State=as.factor(State_abb)) %>%
  ggplot(aes(x=week_num,color=State)) +
  stat_summary(aes(y=group_14_pre,linetype="solid",color=State),geom = 'line') +
  stat_summary(aes(y=group_14_post,linetype="dashed",color=State),geom='line') +
  stat_summary(aes(y=pre_3,linetype="longdash"),geom='line',color="black",lwd=1.4) +
  stat_summary(aes(y=pre_2,linetype="dotted"),geom='line',color="black",lwd=1.4) +
  geom_vline(aes(xintercept = 13),color="grey",linetype="dashed",lwd=1.3)+
  theme_minimal()+
  theme(legend.position = "bottom",legend.title = element_text(color = "black", size = 9),legend.text = element_text(color = "black", size = 7))+
  guides(size = FALSE)+
  xlim(8,16)+
  scale_color_discrete(name="")+
  scale_linetype_discrete(name="",labels=c("Never","Late"),breaks=c("longdash","dotted"))+
  xlab("Week Number")+
  ylab("Mobility")+
  ggtitle("Week 14")


plot4 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_15_pre=ifelse(TreatmentWeek_num2==15&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_15_post=ifelse(TreatmentWeek_num2==15&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(early_late_never==2&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  mutate(State=as.factor(State_abb)) %>%
  ggplot(aes(x=week_num,color=State)) +
  stat_summary(aes(y=group_15_pre,linetype="solid",color=State),geom = 'line') +
  stat_summary(aes(y=group_15_post,linetype="dashed",color=State),geom='line') +
  stat_summary(aes(y=pre_3,linetype="longdash"),geom='line',color="black",lwd=1.4) +
  geom_vline(aes(xintercept = 14),color="grey",linetype="dashed",lwd=1.3)+
  theme_minimal()+
  theme(legend.position = "bottom",legend.title = element_text(color = "black", size = 9),legend.text = element_text(color = "black", size = 7))+
  xlim(8,16)+
  scale_color_discrete(name="")+
  guides(size = FALSE)+
  scale_linetype_discrete(name="",labels=c("Never"),breaks=c("longdash"))+
  xlab("Week Number")+
  ylab("Mobility")+
  ggtitle("Week 15")

library(cowplot)
plot_grid(plot1, plot2,plot3,plot4, labels = NULL,align = "hv")

#ggsave("Figures PoP Appendix/ParallelTrendsTreatmentWeekStates.png")






#================================================================


# Create Parallel Trends plot for Mobility Trneds for All States


data_weekly_series_1 %>% filter(is.na(State_abb)==FALSE)  %>%mutate(cmi_counterf=ifelse(TreatmentWeek_num2==17,cmi_weekly,NA))%>% ungroup() %>%group_by(week_num)%>% mutate(cmi_counterf2=mean(cmi_counterf,na.rm=TRUE))%>% ungroup()%>%mutate(cmi_weekly_pre=ifelse(week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>% mutate(cmi_weekly_post=ifelse(week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%mutate(States=as.factor(State_abb))%>% 
  mutate(TreatmentWeekDummy=as.factor(TreatmentWeek_num))%>%ungroup() %>% mutate(TreatmentWeekLine=ifelse(TreatmentWeek_num2==17,18,TreatmentWeek_num2)) %>%
  ggplot(aes(x=week_num)) +
  stat_summary(aes(y=cmi_weekly_pre,linetype="solid"),geom = 'line',size=0.7,color="red") +
  stat_summary(aes(y=cmi_weekly_post,linetype="dashed"),geom='line',size=0.7,color="red") +
  stat_summary(aes(y=cmi_counterf2,linetype="longdash"),geom='line',size=0.7,color="grey") +
  geom_vline(aes(xintercept = (TreatmentWeekLine-1)),color="black",linetype="dashed")+
  theme_bw()+
  theme(legend.position="bottom",legend.box="vertical",element_blank())+
  labs(
    y = "Mobility", x = "Week Number") + 
  xlim(8,16)+
  scale_linetype_discrete(name="Treatment Status",labels=c("Post","Never","Pre"))+
  facet_wrap(~ States)

#ggsave("Figures PoP Appendix/ParallelTrendsAllStates.png")





#============================================================


# Create Balance Table

data_balance_1 <- data_weekly_series_1 


covs2 <- data_balance_1 %>% ungroup()%>% mutate(`GOP win`=gop_win_2016,`Republican Per.`=gop_per_2016,`Social Capital`=sk2017_alt,Turnout=Turnout12and16,`Republican Governor`=Republican_gov_2020,`Vaccination Per.`=analysis_value,PerCapitaInclog=log(PerCapitaInc),PopDensity2010log=log(PopDensity2010+1),`Unemployment`=UnempRate2018,`Education`=Ed5CollegePlusPct,`Pct White`=WhiteNonHispanicPct2010,`Pct Black`=BlackNonHispanicPct2010,`Pct Asian`=AsianNonHispanicPct2010,`Pct over 65`=Age65AndOlderPct2010,`Rur.- Urb.`=RuralUrbanContinuumCode2013,`Pct Agr.`=PctEmpAgriculture,`Pct Manuf.`=PctEmpManufacturing,`Pct Serv.`=PctEmpServices,`CMI 2019`=cmi_weekly_2019,`Shelter Pct 2019`=sheltered_in_place_2019) %>%
  mutate(`Pop. Dens.`=PopDensity2010log,`Per Capita Inc`=PerCapitaInclog)%>% 
  select(`Republican Per.`,`Social Capital`,Turnout,`Republican Governor`,`Vaccination Per.`,`Per Capita Inc`,`Unemployment`,Education,`Pct over 65`,`Pct White`,`Pct Black`,`Pct Asian`,`Pop. Dens.`,`Rur.- Urb.`,`Pct Agr.`,`Pct Manuf.`,`Pct Serv.`,`CMI 2019`,`Shelter Pct 2019`) 


balTab1 <- bal.tab(covs2,treat = data_balance_1$post_treatment3,s.d.denom = "treated",continuous = "std",binary="std")

love.plot(balTab1,quick=TRUE,stars = "raw",title = "",threshold = 0.1,position = "none")

#ggsave("Figures PoP Appendix/balance_1_final.png")


#================================================================

# Create Goodman-Bacon Decomposition
library(bacondecomp)

bacon <- bacon(cmi_weekly~post_treatment3,data=data_weekly_series_1[is.na(data_weekly_series_1$post_treatment3)==FALSE&is.na(data_weekly_series_1$cmi_weekly)==FALSE&is.na(data_weekly_series_1$FIPS.y)==FALSE&is.na(data_weekly_series_1$week_num)==FALSE,],id_var = "FIPS.y",time_var = "week_num")

ggplot(bacon) +
  aes(x = weight, y = estimate, color = factor(type)) +
  theme_bw()+
  labs(x = "Weight", y = "Estimate", color = "Type") +
  geom_point(size=2)


#ggsave(filename = "Figures PoP Appendix/Bacon2.png")


#================================================================


# Create Mobility Trends by Vaccination Rate Figure


data_weekly_series_1 %>% filter(is.na(vaccination_bin)==FALSE)%>% ggplot(aes(TreatmentWeek_rel2, cmi_weekly,color=as.factor(vaccination_bin))) +
  stat_summary(geom = 'line',size=1.1) +
  geom_vline(xintercept = -1,color="purple",linetype="dashed") +
  theme_minimal()+
  scale_color_discrete(name="Vacc. Rate",labels=c("Low","High"))+
  xlim(-5,4)+
  xlab("Week Number Relative to Treatment")+
  ylab("Mobility")

#ggsave("Figures PoP Appendix/ParallelTrendsbyVaccinationRate.png")





#================================================================



# Create Mobility Trends by Turnout Figure


data_weekly_series_1 %>% filter(is.na(Turnout12and16_bin)==FALSE)%>% ggplot(aes(TreatmentWeek_rel2, cmi_weekly,color=as.factor(Turnout12and16_bin))) +
  stat_summary(geom = 'line',size=1.1) +
  geom_vline(xintercept = -1,color="purple",linetype="dashed") +
  theme_minimal()+
  scale_color_discrete(name="Turnout",labels=c("Low","High"))+
  xlim(-5,4)+
  xlab("Week Number Relative to Treatment")+
  ylab("Mobility")

#ggsave("Figures PoP Appendix/ParallelTrendsbyTurnout.png")




#================================================================

# Create Correlation Figure of Trust in Others (MRP) and Social Capital
MRP_weekly_series_1 %>% mutate(SocCap=as.numeric(StateSocCap))%>% group_by(State_abb) %>% slice(1)%>%
  ggplot(aes(x=SocCap,y=TrustOthers))+
  theme_bw()+
  geom_jitter()+
  geom_smooth()+
  xlab("Social Capital")+
  ylab("Trust Others MRP")



#================================================================

# Create Correlation Figure of Trust in Government (MRP) and Partisanship

MRP_weekly_series_1 %>% mutate(Republican=as.numeric(rep_per))%>% group_by(State_abb) %>% slice(1)%>%
  ggplot(aes(x=Republican,y=NoCorruption_num))+
  theme_bw()+
  geom_jitter()+
  geom_smooth()+
  xlab("Rep. Vote Share")+
  ylab("Trust Gov MRP")+
  xlim(0.25,0.75)





#================================================================



# Create Plot for Main Effects of Interest for Models with Lag and Lead and Lead Treatment Dummies


data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus6=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-46),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus5=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-39),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus4=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-32),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus3=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-25),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus2=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-18),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus1=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-11),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac1=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(3),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac2=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(10),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac3=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(17),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac4=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(24),1,0,missing = 0))

panel_plac_new1 <- felm(cmi_weekly~post_treatment_placminus6+post_treatment_placminus5+post_treatment_placminus4+post_treatment_placminus3+post_treatment_placminus2+post_treatment_placminus1+post_treatment3+post_treatment_plac1+post_treatment_plac2+post_treatment_plac3+post_treatment_plac4|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)



estimates <- as.data.frame(cbind(rep(NA,5),rep(NA,5),rep(NA,5)))


for (i in 1:5){
  estimates[i,1] <- as.vector(panel_plac_new1$coefficients)[(i+4)]
  estimates[i,2] <- as.vector(panel_plac_new1$coefficients)[(i+4)]+panel_plac_new1$cse[(i+4)]
  estimates[i,3] <- as.vector(panel_plac_new1$coefficients)[(i+4)]-panel_plac_new1$cse[(i+4)]
}





colnames(estimates)[1:3] <- c("Point_Estimate","Lower_CI","Upper_CI")

estimates$Week <- c(-2,-1,0,1,2)
estimates$model <- rep(1,5)


panel_plac_new2 <- felm(cmi_weekly~post_treatment_placminus6*gop_per_2016+post_treatment_placminus5*gop_per_2016+post_treatment_placminus4*gop_per_2016+post_treatment_placminus3*gop_per_2016+post_treatment_placminus2*gop_per_2016+post_treatment_placminus1*gop_per_2016+post_treatment3*gop_per_2016+post_treatment_plac1*gop_per_2016+post_treatment_plac2*gop_per_2016+post_treatment_plac3*gop_per_2016+post_treatment_plac4*gop_per_2016|State_abb+week_name|0|State_abb,data = data_weekly_series_1)

estimates2 <- as.data.frame(cbind(rep(NA,5),rep(NA,5),rep(NA,5)))


for (i in 1:5){
  estimates2[i,1] <- as.vector(panel_plac_new2$coefficients)[(i+16)]
  estimates2[i,2] <- as.vector(panel_plac_new2$coefficients)[(i+16)]+panel_plac_new2$cse[(i+16)]
  estimates2[i,3] <- as.vector(panel_plac_new2$coefficients)[(i+16)]-panel_plac_new2$cse[(i+16)]
}





colnames(estimates2)[1:3] <- c("Point_Estimate","Lower_CI","Upper_CI")

estimates2$Week <- c(-2,-1,0,1,2)
estimates2$model <- rep(2,5)



panel_plac_new3 <- felm(cmi_weekly~post_treatment_placminus6*sk2017_alt+post_treatment_placminus5*sk2017_alt+post_treatment_placminus4*sk2017_alt+post_treatment_placminus3*sk2017_alt+post_treatment_placminus2*sk2017_alt+post_treatment_placminus1*sk2017_alt+post_treatment3*sk2017_alt+post_treatment_plac1*sk2017_alt+post_treatment_plac2*sk2017_alt+post_treatment_plac3*sk2017_alt+post_treatment_plac4*sk2017_alt|State_abb+week_name|0|State_abb,data = data_weekly_series_1)

estimates3 <- as.data.frame(cbind(rep(NA,5),rep(NA,5),rep(NA,5)))


for (i in 1:5){
  estimates3[i,1] <- as.vector(panel_plac_new3$coefficients)[(i+16)]
  estimates3[i,2] <- as.vector(panel_plac_new3$coefficients)[(i+16)]+panel_plac_new3$cse[(i+16)]
  estimates3[i,3] <- as.vector(panel_plac_new3$coefficients)[(i+16)]-panel_plac_new3$cse[(i+16)]
}





colnames(estimates3)[1:3] <- c("Point_Estimate","Lower_CI","Upper_CI")

estimates3$Week <- c(-2,-1,0,1,2)
estimates3$model <- rep(3,5)


estimates_combined <- rbind(estimates,estimates2,estimates3)


ggplot(data=estimates_combined, mapping = aes(x = Week, y = Point_Estimate, ymin = Lower_CI, ymax = Upper_CI,color=factor(model))) +
  geom_pointrange(size = 0.8) +
  theme_bw()+
  geom_hline(yintercept = 0, colour="brown4") +
  scale_color_discrete(name="Table",labels=c("Table 2","Table 3", "Table 6"))+
  xlab("Week Relative to Treatment") +
  ylab("Coefficient Magnitude")





#================================================================


# Create Parallel Trends Sensitivity Analysis (follwing Rambachan and Roth 2019)
library(HonestDiD)

data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus6=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-46),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus5=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-39),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus4=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-32),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus3=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-25),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus2=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-18),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_placminus1=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(-11),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac1=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(3),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac2=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(10),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac3=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(17),1,0,missing = 0))
data_weekly_series_1 <- data_weekly_series_1 %>% mutate(post_treatment_plac4=if_else((as.Date(ref_dt)-Stayorder_Date_date)>(24),1,0,missing = 0))

panel_plac_new1 <- felm(cmi_weekly~post_treatment_placminus6+post_treatment_placminus5+post_treatment_placminus4+post_treatment_placminus3+post_treatment_placminus2+post_treatment_placminus1+post_treatment3+post_treatment_plac1+post_treatment_plac2+post_treatment_plac3+post_treatment_plac4|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)

numPostPeriods <- 5
numPrePeriods <- 6

l_vec <- basisVector(1,5)


DeltaSD_RobustResults <- createSensitivityResults(betahat = panel_plac_new1$beta,sigma=panel_plac_new1$clustervcv,numPrePeriods = numPrePeriods,numPostPeriods = numPostPeriods,l_vec = l_vec)

originalResults <- constructOriginalCS(betahat = panel_plac_new1$beta,sigma=panel_plac_new1$clustervcv,numPrePeriods = numPrePeriods,numPostPeriods = numPostPeriods,l_vec = l_vec)


createSensitivityPlot(robustResults = DeltaSD_RobustResults,originalResults = originalResults)

#ggsave("Figures PoP Appendix/sensitivity_lagleads.png")





#==========================================================================================================================
#======================================================================================================================